%alpha=1;
qx=1./(1-etaI*ixx);

for j=2:flagstar
dqx(j)=(qx(j)-qx(j-1))/ds;
end

for j=1:flagstar
s=ss(j);
betas=beta0+beta1*s; 
lambdas=max(lambda0*(1-zeta0*s^zeta1),0);
zetas=max(eta0*(1-s^eta1),0);
zetan(j)=zetas;

ix=ixx(j);
x=xx(j);
phii=ix-etaI*ix^2/2-deltaI;
gx(j)=phii-lambdas/(betas+1)-zetan(j)*(qx(j)-qxB(j))/qx(j);
cx(j)=A-ix-x;

rp=gamma*sigma^2+lambdas*(betas/(betas-gamma)-betas/(betas+1-gamma))-lambdas/(betas+1)+zetan(j)*(qx(j)-qxB(j))/qx(j)*((qxB(j)/qx(j))^(-gamma)-1);
phiX(j)=xx(j)/ss(j)-etaX/2*xx(j)^2/ss(j)^2-deltaX;
phiI(j)=ixx(j)-etaI/2*ixx(j)^2-deltaI;

r=cx(j)/qx(j)+gx(j)-rp+(phiX(j)-phiI(j))*ss(j)*dqx(j)/qx(j);
mx(j)=xx(j)/alpha;
thetaU=x*(1-etaI*ix);
thetaS=-(1-alpha)*mx(j)/qx(j);

qU(j)=(A-ix+(phiX(j)-phiI(j))*ss(j)*dqx(j))/(r+rp+thetaU-phii+lambdas/(betas+1)+zetan(j)*(qx(j)-qxB(j))/qx(j));
qS(j)=(A-ix-x+(phiX(j)-phiI(j))*ss(j)*dqx(j))/(r+rp-phii+lambdas/(betas+1)+zetan(j)*(qx(j)-qxB(j))/qx(j));
rpx(j)=rp;
rm(j)=r+rp;
rx(j)=r;

rpU(j)=rp+thetaU;
rpS(j)=rp+thetaS;

%ux(j)=((r+rp-phii+lambda/(betas+1))*rho^(-psi))^(1/(1-psi));
ux(j)=((cx(j)/qx(j))*rho^(-psi))^(1/(1-psi));
end